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Abstract 

Time-delayed control in a balancing problem may be a nonsmooth function for 
a variety of reasons. In this paper we study a simple model of the control of an 
inverted pendulum by either a connected movable cart or an applied torque for 
which the control is turned off when the pendulum is located within certain regions 
of phase space. Without applying a small angle approximation for deviations 
about the vertical position, we see structurally stable periodic orbits which may 
be attracting or repelling. Due to the nonsmooth nature of the control, these 
periodic orbits are born in various discontinuity-induced bifurcations. Also we 
show that a coincidence of switching events can produce complicated periodic and 
aperiodic solutions. 
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c3 ■ 1 Introduction 

The subject of balance has received considerable recent attention. Local measurements 
of muscles [U El E] and experiments that highlight the influence of vision pEJ |5] have 
led to an improved understanding of the key physiological elements in human balancing 
tasks. From a theoretical perspective, progress has been made in analyzing systems 
with time-delay [6j. Time-delay is used to model the reaction time of the controlling 
mechanism and is a near ubiquitous element of mathematical models of balancing tasks. 
A current challenge is to incorporate new experimental observations into mathematical 
models and interpret the results. Time-delayed balance control is also a fundamental 
problem in robotics [7J EJ |9j, however the control strategies used in engineering are 
typically distinct from those identified in physiology [ID] . 
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Time-delayed PD control (for which the applied forces are determined from mea- 
surements of position (P) and its derivative (D)) has been used in models of the control 
of a stiff beam or rod that is attached at its base to a controllable cart [TTJ [T2] and 
stick balancing [T5| . In both cases for small time-delay it is possible to choose control 
parameters so that the model is successfully directed to the vertical position. With an 
increase in the delay time, for fixed control parameters, the vertical position becomes 
unstable which may result in stable oscillations about the vertical position. Beyond a 
critical value of the delay, the control is incapable of stabilizing the system at the vertical 
position. A similar response has been found with other smooth control laws [2]. Some 
authors have used more complex models to incorporate additional physical features such 
as friction [15]. In general, time-delayed control models are inherently difficult to ana- 
lyze because they are infinite-dimensional. However, dynamics near local bifurcations 
of delay differential equations are well described by low-dimensional systems of ordinary 
differential equations derived via a centre manifold analysis [Ml IT7] . 

Recently several researchers have proposed nonsmooth balancing models. Exper- 
imental observations of human balance in quiet standing [18] and studies of human 
balancing tasks [191 [20] suggest a switching control due to intermittent muscle move- 
ments [21] . Specifically, experiments suggest a state-dependent control such as the "drift 
and act" method of fTS] which simply turns the control off when the system nears equi- 
librium. Muscle control may be active or passive, the latter refers to muscles which 
intrinsically resist motion away from equilibrium. The presence of control is potentially 
detrimental to obtaining equilibrium if the system is naturally approaching equilibrium; 
in [22] the authors consider a switching control that lessens this effect. In [23], the 
authors present a control that acts only after waiting for a time longer than the delay in 
order to gain sufficient information from the system to be able to perform more effective 
control. Hysteretic control laws, which are common in temperature control, have been 
considered in balancing models [2~i|l2"5]. In mechanical systems small gaps between gears 
create backlash which another source of nonsmoothness and in the absence of friction 
eliminates the possibility of perfect stable equilibrium at the vertical position [26]. From 
an engineering or robotics viewpoint, a switching control may require less cost and be 
easier to implement mechanically or may be necessary due to a nonzero sampling time. 
An additional benefit is that nonsmooth control laws are often able to stabilize the sys- 
tem for arbitrarily large delay in simple mathematical models. For an introduction to 
switching in control systems, see for instance [27] . 

Mathematical models that incorporate both time-delay and switching conditions are 
usually particularly difficult to analyze, yet there are common mechanisms that induce a 
transition from simple to complex dynamics and are characteristic of such systems. For 
instance a periodic orbit may develop a tangency with a switching condition, or undergo 
switching at times that differ by exactly the delay time of the system [2l| |28j [29]. Both 
scenarios correspond to a co dimension-one bifurcation of the periodic orbit. 

The purpose of this paper is to investigate the effect of the application of control 
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with two different switching rules on simple balancing models. We study nonlinear 
equations of motion for a pendulum combined with time-delayed PD control by a force 
applied either by a movable cart or as a torque. We consider two different ON/OFF 
switching rules for the application of the control. For the majority of this paper we 
analyze a switching rule that involves both position and velocity, a later section of the 
paper is devoted to a switching rule based on position only. The goal is to reveal and 
understand novel dynamics resulting from the combined effects of time-delay, switching 
and nonlinearity which are generic and therefore expected to be prevalent in a wide 
range of balancing systems. 

The remainder of this paper is organized as follows. The equations of motion inves- 
tigated are stated in §2.11 In §2.21 we introduce two switching rules that divide phase 
space into various ON and OFF regions. We then focus on the first switching rule. 
Section 12.31 summarizes dynamics when the delay time is zero. In this case orbits may 
become "stuck" to a switching manifold. In the presence of small time-delay this sliding 
motion becomes rapid switching motion, which we refer to as zigzag motion about the 
switching manifold, §3.11 This motion corresponds to a jittery motion of the pendulum 
restricted to one side of the vertical position maintained by an intermittent application 
of the control. The time-delay may alternatively induce spiral motion corresponding to 
oscillations about the vertical position, but for short time-delay zigzag dynamics dom- 
inate §3.21 The bifurcation structure is detailed in §3.31 In particular we prove that 
both stable and unstable zigzag periodic orbits may be born in discontinuity-induced 
bifurcations. These bifurcations are described formally through asymptotic expansions 
in §3.41 Homoclinic zigzag orbits are the subject of §3.51 

In §H we consider longer values of time-delay for which numerics reveal complex 
dynamics that are not explained by the small delay asymptotic expansions. We describe 
a novel bursting-like attractor that exhibits different behaviours on two distinct time- 
scales, §4.11 In §4.21 we map out the stability of the pendulum at the vertical position 
in the plane of the control parameters by linearizing the equations of motion. 

The switching rule based on position only is studied in $5] This rule corresponds 
to the controller neglecting control when near the vertical position. We again com- 
pute asymptotic expansions and prove the existence of stable periodic orbits. Finally a 
summarizing discussion is presented in §|6j 

2 Simple Balancing Models 

Here we detail the mathematical models studied in this paper. To give our results wide 
applicability we have chosen to study dimensionless, inverted pendulum-based equations 
of motion that have been used as models in both human balancing tasks and mechanical 
systems. For simplicity we consider all motions to be restricted to a plane and ignore 
both friction and noise. Section I2TT1 introduces the equations of motion and PD control, 
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§2.21 details two ON/OFF control mechanisms and describes basic dynamics. Dynamics 
in the absence of delay are described in §2.31 



2.1 Equations of Motion 

The task of vertically balancing a long stick has been modelled as an inverted pendu- 
lum with an applied torque [10J [231 [13] . Upon an appropriate time scaling and other 
reductions the equation of motion may be written simply as 

6 - sin(0) = F , (1) 

where is a dimensionless quantity representing the angular displacement of the stick 
from vertical and F denotes pivot control due to, say, the finger or hand of the human 
performing the balancing actions. 

Equation flTJ also provides a simple model of human postural sway, where F repre- 
sents ankle torque [22] . For postural sway it is important to treat non-control compo- 
nents of F. In particular, ankle torque has an intrinsic passive stiffness which provides 
some stability but is regarded as inadequate to maintain quiet standing [H [2] . The ankle 
joint also provides damping. If human postural sway is overdamped it may be suitable 
to model the motion with a first order differential equation P3H EH]- Analyses of (jTJ 
provide a basis for more complex motions such as balancing with hip movements and 
bipedal walking in humans and robots. 

Planar dynamics of a vertical rod controlled by a moving cart have been modelled 

by 

1 _ ¥H COS 2 (0)^ + — 2 sin(20) - sin(0) = F cos(0) , (2) 
4 / 8 

[TTl [T2"l . where m = M pcnduium — denotes the fraction of the mass of the system that 

1 ' " -Mpendulum+Mcart J 

belongs to the pendulum and friction is ignored. If the cart is much more massive than 
the pendulum, i.e. M cart M pendulum , then m = is a useful approximation and the 
equation of motion simplifies to 

0-sin(0) = Fcos(0) . (3) 

Equation (T5]) exhibits dynamics similar to ([3]) for small values of m. A cart model with 
friction is studied in [15]. 

In this paper we study both (CQ) and (j3J). For convenience we let <fi = 6 and write 

= 0, 

• (4) 
(f) = sin(0) + FG(6) , 

where G{9) = 1 or G{9) = cos(0). For physical scenarios that involve only small changes 
in angular displacement, it is suitable to linearize equation of motions in 0. We do not 



4 



use a linear (small angle) approximation in 8 in order to investigate invariant solutions 
created in bifurcations at 8 = 0. 

The control force, F, is a time-delayed function of 6 and (p. In the context of human 
balancing tasks time-delay models neural transmission time. One of the simplest control 
laws is PD control [HUE]: 

F N = a6(t - r) + b<f>{t - r) , (5) 

where a and b are scalar control parameters and r > is the delay time. For inverted 
pendulum-type problems with PD control that is continuously applied, the vertical po- 
sition is typically stable for control parameters that lie in a roughly D-shaped region in 
the (a, 6)-plane [T21T13J. Pitchfork bifurcations and Hopf bifurcations form the boundary 
of this region. Outside the region there may exist stable oscillations about the vertical 
position, complex dynamics or even failure for the system to attain a physically mean- 
ingful bounded solution. The area of the D-shaped region reduces as the value of r is 
increased. In [17], the authors study fl2]) with <^ and find that the D-shaped region 
shrinks to a point at some critical value of r. 



2.2 ON/OFF control 

By ON/OFF control we mean simply that at some times the control is implemented, 
whereas at other times the control is absent, i.e. F — 0. In this paper we study two dif- 
ferent ON/OFF controls based on position in the (8, (^)-plane. In contrast the ON/OFF 
mechanism in the so-called "act-and-wait" strategy [231 E2J S3] is based on the time 
elapsed. When control is applied we use the PD control (jSJ). We refer to (]3]) with F = 
as the OFF system and (]3]) with (j3J) as the ON system. 
The two control laws we analyze are 

Fon , 8(t - t) (<f>{t - t) - s8{t - r)) > 
, otherwise 

i^ON, |0(*-T)|><7 (7) 

, otherwise 

where s < and o > 0. Equations (jSJ) and flTJ) are respectively taken from [22] and 
[T§1 El] , which focus on human postural sway. The control law (jB]) defines two switching 
manifolds 

Vi = {(e,se)\eeR}, 

5 2 = {(O,0) U 

that divide the (8, </>)-plane into four regions. We refer to these regions as ON and OFF 
regions, as indicated in Fig. [TJ Similarly (J7J) defines two switching manifolds 

5 3 = {( C r,0) | 0eR}, 

£ 4 = {(-<r,0) | U 
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that divide the (9, <^)-plane into an OFF region that is a vertical strip centred about 
9 = 0, and two ON regions, Fig. [TUJ The system dl])-© retains the usual symmetry by 
(9, 4>) i — y —(9, (f>) with either switching condition. 

The system flU)-© with either flS} or (J7]) is a piecewise-smooth, discontinuous, delay 
differential equation system. Due to the presence of time-delay the phase space of 
the system is infinite-dimensional [61 [351 EH], nevertheless it is convenient to picture 
dynamics in the (9, 0)-plane. It is difficult to consider all possible initial conditions 
(which are curves, (9(t),(f)(t)) for t G [— r, 0]). However, when F = 0, (J1J reduces to a 
two-dimensional ODE system (the OFF system). Consequently, if the point (9(0), (f)(0)) 
lies in an OFF region and the trajectory (9(t),<p(t)) has been in this OFF region for 
a time equal to at least r, then F = for all t G [0, r] and so the trajectory at any 
positive time is independent of its location at any time prior to t = 0. Hence, with 
this assumption, the initial condition may be thought of as merely the location of the 
trajectory at t = 0. Since it is straight-forward to understand the dynamics of the OFF 
system, we consider as initial conditions only points on switching manifolds at which 
the vector field of the OFF system points into the neighbouring ON region. Admittedly 
this restriction omits some behaviour of the system, but we believe it captures all the 
important and physically meaningful dynamics. 




Figure 1: The (9, 0)-plane for (S])-(j5]) with (|6]). Two trajectories were numerically 
computed when s = —0.3, r = 0.5, (a,b) = (1.5,4) and G(9) = cos(0). One trajectory 
zigzags about Si and tends to the origin, the other trajectory spirals out from the origin 
limiting upon a stable periodic orbit. W s is the stable manifold of the origin for the 
OFF system. 
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2.3 Dynamics in the Absence of Delay 



Here we describe the behaviour of (j4])-(l6]) when r = 0. We begin with brief analyses of 
the individual OFF and ON systems when r = 0, then describe the effects of switching 
rule ([6]). 

The OFF system 

The OFF system, given by (j3J) with F — 0, represents a classical inverted pendulum 
and is independent of the time-delay. The system is Hamiltonian and the function 



is a suitable Hamiltonian function for this system (i.e. H(8(t), 4>{t)) is constant for any 
solution). In §3 .41 we will use (flOl) to measure the variation of trajectories in the presence 
of small delay over the course of individual zigzag oscillations. The equilibria of the OFF 
system are (mr, 0), for n G Z. If n is even, the equilibrium is a saddle; if n is odd, the 
equilibrium is a centre. 

The ON system 

In the absence of delay, the ON system for (jl])-© is 



which corresponds to instantaneous PD control. For both choices of G, the origin is 
an equilibrium of (12) and it can be characterized with a standard stability calculation. 
For a < 1 the origin is a saddle, otherwise it is a node or a focus. The node or focus is 
stable [unstable] for b > [b < 0]. 

For G(9) = cos(#), (TlTl) has infinitely many equilibria, and for a > 1 the equilibrium 
(±#* os ,0) is a saddle, where 9* os denotes the smallest positive 9- value for the equilibria. 
When b > 0, the value a = 1 corresponds to a subcritical pitchfork bifurcation at the 
origin. 

For G(9) = 1, the origin is the only equilibrium of fTTTl) for a > 1. For < a < 1 
the equilibrium (±#*,0) is stable when b > 0, where 9\ denotes the smallest positive 
#-value for the equilibria. Consequently, when b > 0, the value a = 1 corresponds to 
a supercritical pitchfork bifurcation at the origin. For both choices of G the pitchfork 
bifurcation forms the left boundary of the D-shaped stability region, §2.11 



H(9,<p) = ^ 2 + C os(6) , 



(10) 



9 = <f>, 

<j> = sin(0) - (aO + b<f>)G(9) , 



(11) 
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The full system 

Here we analyze (H])-© when r = 0, which may be written as 



sin(0) - {a9 + hj>)G{9) 

<t> 

sin(0) 



- s9) > 

(12) 

otherwise 



This system is a Filippov system [37J EE] because it is discontinuous on the switching 
manifolds, Si and S 2 ©• Dynamical behaviour that lies entirely within either an ON 
or OFF region is determined purely by either the ON or OFF system. Thus here we 
focus on trajectories that impact a switching manifold. For impacts on S 2 , since 6 = <f> 
for both ON and OFF systems, trajectories simply arrive at S 2 from the neighbouring 
OFF region and immediately enter the adjacent ON region. 

For impacts on Si, we identify regions where the trajectory remains on Si for some 
time. A section of Si along which the ON vector field points into the OFF region and 
the OFF vector field points into the ON region is known as an attracting sliding region 
[3SJ EH]- A trajectory that arrives at an attracting sliding region becomes stuck on the 
switching manifold and slides. The manner by which sliding dynamics evolve is usually 
defined by Filippov's method, [37J [381 EHJ SO], which we explain below. 

First let us locate sliding regions on Si. At an end of a sliding region the vector field 
of either the OFF system or the ON system is tangent to Si. Such a point is referred 
to as a grazing point and a trajectory that exhibits this tangency is a grazing trajectory. 
The slope of the vector field of the OFF system on Si, | = ^f> is tangent to Si when 

sin(#) 9 

m = 9 ~ S ° • (13) 
Similarly, the slope of the vector field of the ON system on Si is tangent to Si when 

g{9) = - s 2 - (a + bs)G{9) = . (14) 

9 

Roots of J 7 and Q correspond to possible grazing points and boundaries of sliding regions. 
For any s G (—1,0], T has a unique root, #2.^J G (0, 7r], and the OFF vector field points 
into the ON region on Si for < 9 < 0g^J. However this point does not influence 
physically meaningful dynamics unless the value of s is close to —1 because whenever 

s>'\[l^ -0.7979, # g ° r L F >|. 

To identify roots of Q, ( fT4l) . we consider the two cases of G(9) separately and omit 
some messy but elementary calculations. For G(9) = 1, there is a unique root for 
^<a + bs + s 2 <l given by 6>°^ z G (0, |). It follows that the subset of Si for which 



8 



gON 



graz 



< 9 < | is an attracting sliding region. For G(9) = cos(#), when a > 1 — bs — s 1 



sliding occurs on Si for < 6 < 6^ z where 9^ z is the smallest positive root of Q. 

Dynamics on an attracting sliding region are governed by the unique convex combi- 
nation of the ON and OFF vector fields that is tangent to the region at each point. We 
write 



where 







" 9 ' 


= (1- 






J. 


slide 


" 9 ' 


and 


' 9 ' 


refer to 


J. 


OFF 


.i. 


ON 



" 9 ' 




" 9 ' 


J. 


+ q 


J. 




OFF 





(15) 



ON 



with F = and ffTTj) respectively, and q is 



(a+bs)6 cos(6) 



slide = S^slide- Upon 

which leads to the 



a ^-dependent scalar quantity determined by the requirement 
substituting (p = s9 into f|T5|) . </> s iide = s^siide yields q — 
explicit solution 

9 s hde(t) _ g 1 
0slide(£) J ° [ S 

Therefore when s < attracting sliding trajectories approach the origin; when s — 0, 
attracting sliding regions are intervals of equilibria. 

Fig. |2]-A illustrates typical dynamics of ( I12p when G(9) = cos(9). There are two 
bifurcations. At a = 1, two saddle equilibria are created that exist to the right of 
this line. This bifurcation is a pitchfork bifurcation of the ON system but not a bona 
fide pitchfork bifurcation of the full system because the origin is a non-differentiable 
point. Throughout the paper we refer to a = 1 as a pitchfork-like bifurcation. At 
a = 1 — bs — s 2 (labeled SL in the figure), a sliding region is created on Si that grows 
in size with increasing a. Panel B summarizes dynamics when G{&) = 1. 



3 Dynamics with Delay 
3.1 Zigzag and Spiral Dynamics 

As discussed in §2.2[ we consider the forward evolution of a point on Si or S2 in 
order to identify the basic behavior of (SJl- (EJ) • Typically the corresponding trajectory 
immediately enters the neighbouring ON region, and by symmetry we may assume that 
this is the ON region with 9 > 0. Notice 9 > whenever <fi > 0, thus the trajectory 
cannot exit the ON region through S 2 . 

One possibility is that the trajectory enters the ON region and remains in the ON 
region for all time. In this case a physically meaningful stable solution is not attained. 
Roughly speaking this occurs for small values of the control parameters. Alternatively 
the trajectory intersects Si for the first time at some t\ > 0. In this case generically the 
trajectory then resides in the OFF region, with 9 > 0, for some nonzero time. Assuming 
the trajectory intersects either Si or S 2 at a later time, let £2 denote the earliest such 
intersection time. 
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If t2 > t\ + r, that is, if the trajectory has been in the OFF region for a time greater 
than or equal to the delay time, r, then the fate of the trajectory at times later than £ 2 
is independent of its location at any time prior to t%. Therefore in this case the location 
of the trajectory aX t = t 2 may be thought of as a new initial point, §2.21 If instead 



A 




0.5 1 1.5 

a 



Figure 2: Bifurcation diagrams of the system (jlJ)-(jEJ) in the absence of delay, ( JT2l) . with 
s = —0.3 and 6 = 2 for G(9) = cos((9) in panel A and G(9) = 1 in panel B. Solid and 
dashed curves denote stable and unstable equilibria, respectively. The dotted curves 
indicate 9™ z , which is a root of G{9) (JHJ) and corresponds to the onset of sliding. PF 
- pitchfork-like bifurcation; SL - grazing point intersects origin. Included are represen- 
tative phase portraits with dash-dot lines indicating the corresponding values of a. The 
system exhibits dynamics similar to that shown here for a significant range of s and b 
values; some bounds on these values are given in the text. 
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t 2 < ti+T, the behaviour of the trajectory is dependent on the control for t e (t 2 , ti+r). 
Trajectories with this property require more effort to analyze beyond t = t\ + r and 
may be particularly complicated but occur commonly only when r is relatively large, 
see mi 

Here we consider the former case, t 2 > t\ + r, and classify two basic types of delay- 
induced dynamics, as noted qualitatively in [22J. At the time t — t\ + r the applied 
control is switched off. It is instructive to consider the location of the trajectory at this 
time, i.e. the point (9{t\ + r), <p(ti + r)), in relation to the stable manifold of the origin 
for the OFF system, W s , shown in Fig. [TJ If the point {9(t\ + r), <p(ti + r)) lies above 
W s then the point (0(t 2 ), 0(t 2 )) lies on S x ; if (0(t a + r),0(ti + r)) lies below W s then 
(6*(^2) , 0(^2)) lies on S 2 . (In the special case that the point (9(ti + r), 0(£i + t)) lies on W s 
(H(6(ti + r), (f){ti + r)) = 1), the trajectory coincides with W s after t\ and never exits 
the OFF region.) If (0(0), 0(0)) and {9{t 2 ),(j)(t 2 )) both lie on S x we refer to the part of 
the trajectory between these two points as a zigzag oscillation. Similarly if {9(0), 0(0)) 
and {9(t 2 ),(p(t 2 )) both lie on S 2 we refer to the same part of the trajectory as half a 
spiral oscillation. Zigzag oscillations often come in succession, as do spiral oscillations. 
We have not observed sustained switching between zigzag and spiral motion. Zigzag 
motion is typical for small values of r and spiral motion is typical for large values of r, 
however, for a carefully tuned combination of the control parameters zigzag and spiral 
motion may coexist, as in Fig. [U 

3.2 Domination of zigzag trajectories for small delay 

Let us consider the forward orbit of a point on Si, (9 ,s9 ), as shown in Fig. [31 for 
the system (HJ-OH]). With a small enough time-delay (r < — s is sufficient), the orbit, 
as governed by the OFF system, does not cross the 0-axis before the control is applied. 
If the applied control is sufficiently large so that 0(£) < in (@|, the orbit abruptly 
changes heading and soon reintersects Si. After a time r beyond this reintersection, 
the control is switched off. As long as the applied control is not so strong that the orbit 
has overshot W s , the orbit then continues back to Si and the process repeats. For all 
time 9(t) is strictly decreasing because 0(t) is always negative. In other words the orbit 
zigzags into the origin. 

Numerical investigations suggest that as long as \ s\ is not too large, say — 0.4 < s < 0, 
and r < — s, then there exists a range of parameters for which the forward orbit of any 
point (6*o, s9 ), with \9 \ < 9* os for G(8) = cos(6>) and \9 \ < | for G(9) = 1, zigzags 
into the origin as in Fig. |3j With stronger control (specifically larger values of a) orbits 
may cross W s producing spiral dynamics. With larger delay the orbits enter the first 
quadrant of the (9, 0)-plane. In this case orbits may still zigzag, but not necessarily 
approach the origin. The next section investigates these dynamics. 



11 



Figure 3: The forward orbit of a point on Ei for (jlj)-(jBD with s = —0.3, r = 0.25 and 
{a,b) = (2.5,2). 



3.3 Bifurcation sets 

In the previous section we argued that if r < — s, then for appropriately chosen a and b 
orbits simply zigzag into the origin. For the remainder of this section we consider small 
values of s so that the condition r < — s may not be satisfied. 

Fig. 0]-A is a numerically computed bifurcation set of dlJ-Q when G{9) = cos(6>). 
For this figure we have fixed the value of s at —0.01; by setting the vertical axis to - 1 - 
the picture is roughly unchanged for different small s < 0. We fixed b = 2 for Fig. H] and 
numerically have observed that the bifurcation structure is qualitatively the same for 
different values of b > 0. The figure shows dynamics only for 9 > 0; identical dynamics 
occurs for 9 < since the system is symmetric. 

For sufficiently small r, and a > 1, Fig. 0]-A reflects the results of §3.21 orbits zigzag 
about Ei and approach the origin. An increase in r leads to the creation of a pair 
of zigzag periodic orbits in a classical saddle-node bifurcation. One periodic orbit is 
stable, the other is unstable. If a > 1.54, a further increase in r destroys the stable 
zigzag periodic orbit via a homoclinic connection to the saddle equilibrium, (9* os ,0). 
Otherwise an increase in r destroys the unstable zigzag periodic orbit by a collision of 
the periodic orbit with the origin. This is a discontinuity-induced bifurcation (labelled 
sub DIB in Fig. HJA) that in some ways resembles a subcritical Hopf bifurcation. Indeed 
the amplitude of the periodic orbit grows at a rate proportional to the square root of 
parameter change (shown by Fig. H]-B) but it does not correspond to the occurrence of 
purely imaginary stability multipliers, nor does the periodic orbit encircle the equilib- 
rium. Furthermore, the system has identical dynamics for 9 < 0, so a second unstable 
zigzag periodic orbit is created simultaneously and exists for 9 < 0. 

The curve of discontinuity-induced bifurcations and the curve of saddle-node bifur- 
cations are tangent to one another at their point of intersection. The criticality of the 
discontinuity-induced bifurcation changes here. That is, along the discontinuity-induced 
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bifurcation curve for a < 1.09, a stable zigzag periodic orbit existing for 9 > is cre- 
ated. Beyond — - = 4 we were unable to numerically continue the bifurcation curve 




Figure 4: Panel A is a bifurcation set of (SJ)-(J6]) when s = —0.01, 6 = 2 and 
G{6) = cos(6>). The solid curves are the result of numerical computations. The dashed 
curves correspond to equations derived in §3.41 and §3.51 PF - pitchfork-like bifurca- 
tion; SN - saddle-node bifurcation of periodic orbits; HC - homoclinic bifurcation; DIB - 
discontinuity-induced bifurcation (super and sub are abbreviations for supercritical and 
subcritical, respectively). Included are six representative sketches of trajectories in the 
(8, </>)-plane. The sketches are exaggerated for clarity - in reality 6{t) changes by an 
extremely small amount over each zigzag. Panel B is a bifurcation diagram correspond- 
ing to the horizontal dash-dot line in panel A (i.e. for r = 0.035). The solid [dashed] 
curves denote stable [unstable] equilibria. The double curves correspond to the maxi- 
mum ^-values of periodic orbits with line style indicating stability in the same fashion. 
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when s = —0.01. Certainly as r increases the discontinuity-induced bifurcation ap- 
proaches the pitchfork-like bifurcation (a = 1). A more detailed investigation of this 
area of parameter space is presented in §4.11 for s = —0.1 and suggests that complex 
dynamics may occur here. 

Numerical simulations indicate that for different values of b > 0, the system exhibits 
a basic bifurcation structure identical to that shown in Fig. |H As b — > + the intersection 
of the HB and SN curves appears to approach (a, — -) = (1, 2) and the HC curve seems 
to limit on the SN curve and a — 1. For the values of b and s corresponding to Fig. HI 
when r = 0, sliding occurs for a > 1 — bs — s 2 = 1.0199 which is so close to a = 1 that 
we have chosen not to indicate it in the figure. 

A bifurcation set for the same parameter values as Fig. H] but with G(9) = 1 is 




Figure 5: A bifurcation set and bifurcation diagram of (SJ)- (JHJ) when s = —0.01, b = 2 
and G(9) = 1. The meanings of the abbreviations and the significance of the line styles 
are the same as in Fig. HI 
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shown in Fig.[5]A. Notice the discontinuity-induced bifurcation curve is unchanged. This 
is because the difference in the two functions of G(9) is 0(8 2 ) and the discontinuity- 
induced bifurcation is local to the origin. The quadratic terms affect the criticality of the 
bifurcation; indeed changing the function G has reversed the criticality on the bifurcation 
curve, Fig. EJ As above, here a curve of saddle-node bifurcations of zigzag periodic orbits 
emanates from the codimension-two point at which the criticality changes. The saddle- 
node bifurcation curve is shown in Fig. [5]A up until the periodic orbit at the bifurcation 
is no longer physically meaningful (where it includes values greater than |). For 
G(9) = 1 the equilibria created at a = 1 are stable and for this reason no homoclinic 
connection forms analogous to that for G(8) = cos(#). 

3.4 Series expansions in s and r 

In this section we derive asymptotic expressions for bifurcations relating to zigzag dy- 
namics that were identified numerically in the previous section in order to gain a greater 
understanding of the bifurcations. The methodology we employ requires the period of 
zigzag oscillations to be small. For this reason we assume that both r and s are small, 
which is consistent with observations that zigzag dynamics occurs for small r and s, and 
obtain explicit expressions for zigzag orbits as series expansions in s, r and t. 

Let T be the forward orbit of a point (9q, s9 ) at t = 0, see Fig. El Let (6>i, 0i) denote 
the location of V at t — r. This is the first switching point of T. Let T int be the next 
time at which T intersects Si, if such an intersection exists. Then the second switching 
point occurs at t = T int + r and we denote this point by (6* 2 , 4>2)- For small s and r it is 
reasonable to assume that the second switching point lies in the OFF region above W s , 
so then T will exit the OFF region through Si at some point (#3, 5^3). Naturally we are 
interested in the difference between 63 and 9q as this indicates whether T is approaching 
or moving away from the origin. Algebraically it is easier to instead compute the change 
in the Hamiltonian, H(9, (ft) ffTQj) . between the two points. Time evolution of the OFF 




Figure 6: A sketch of the forward orbit, T, of a point, (#oi s9q), on Si. 
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system does not change the Hamiltonian, so it is equivalent to look at 

AH = H 2 -H 1 , (16) 

where 

H 1 = H(6 1 ,<f>i), H 2 = H(9 2 ^ 2 ). 

We now compute the first few terms of AH as a series expansion in s and r. To 
derive the expansion we first directly use the governing differential equations (J4j) — (J6J) to 
express T as a series in s, r and t, then solve for Ti nt , and finally evaluate Hi and H 2 . 
Since the system under consideration includes time-delayed switching, for our purposes 
this approach is preferable to expanding r within the differential equations with the idea 
of reducing the delay differential equations to a r-dependent ODE system. A note with 
regards to notation: we use 0(\ ■ \ k ) to denote terms of order k or higher in the given 
variables. We assume that s and the period of oscillations are both 0(r). 

The orbit V is the solution to the initial value problem, (SJ)- (JHJ) with {0(0), 0(0)) = 
(9o,s9q). For t G [0,r], T is governed by the OFF system. By substituting a series of 
the form 9(t) = Y2i Ylj c ij{0o)s l P into (j3J) with F = and solving for the coefficients we 
obtain 

9(t) = 9 + o st + l - sin(# )t 2 + 0(\s, t| 4 ) , (17) 
which is valid for t G [0, r\. Substituting t = t into (ITTj) yields 

9 X = + 9 sr+^sm(0 )T 2 + O(\s,T\^ , (18) 
0x = 9 s + sm(9 )T + O(\s,T\ 3 ) . (19) 

Proceeding in a similar fashion we obtain an explicit expression of 9(t) for t > r. We 
provide further comments on the derivation below. The expansions are centred about 
(6*i, (pi) instead of (9 , s9 ) because this provides some simplification and then powers of 
t — t naturally appear. For small s, t and t and assuming T int = O(r), T is given by 

' 0i + [a0 x + sin(^)T) (t - t) + | sin(0i)(i - r) 2 + 0(\s, r, t| 4 ) , t G [0, r] 

9(t ) = J 01 + ( s0 i + sin(0i)r) (t - r) + (a 3 + a 4 s) (t - r) 2 + a 6 (t - r) 3 + 0(|s, r, t| 4 ) , t G [r, 2r] 

U j 9x + aiT 3 + (s0i + sin(6»i)r + a 2 r 2 ) (t - r) + (a 3 + a 4 s + a 5 r) (t - r) 2 + a 6 (i - r) 3 + 0(\s, r, i 

tG[2r,T int + 
(20) 
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where 



&A0) = — ab9G(9? , 
6 

a 2 {9) = ^ab9G(9) 2 , 



a z {9) = --(a0G(0)-sin(0)) , 



a 4 (9) = ~b9G(9), 
a 5 {9) = -^ab9G(8) 2 , 



a 6 (9) = -hsin(9)G{9) 2 , 
o 

a 6 (9) = -bG{9){a9G{9) - sin 



(0)) 



and in (T2"Uj) the a's are evaluated at 9 = 9\. 

The top-most expression of fl20l) follows from combining ffl7]) and ffl8|) . The middle 
expression of (120]) is obtained from the ON system and substituting the top-most ex- 
pression in place of 9(t — r) and using the time derivative of the top- most expression for 
(f){t — t). For this reason the middle expression is valid only for t G [r, 2r]. Then using 
the middle expression for the time-delayed components, 9(t — r) and <p(t — r), we obtain 
the bottom expression which is valid for t G [2r, 3r] , if and only if T; nt > 2r. Continuing 
in this fashion one may build up an explicit expression for T to any desired order in s, 
r and t with formulae that are valid in [nr, {n + l)r], with n G Z, until i = T int + r, 
beyond which the OFF system governs V for some time. The solution to the ON system, 
9(t), is one degree more differentiable than 9(t — r). Consequently, 9(t) is C n_1 at each 
t = nr with n£Z and nr < T; nt + r. Therefore the series solution to V for t G [3r, At] 
differs from the bottom expression of (T20T) only in terms that are 0(\s, r, £| 4 ). Since ( |20|) 
has no explicitly stated quartic or higher order terms the bottom expression is valid for 
all t G [2t, T int + t] as stated. 

The intersection time, T int , is defined by 0(T int ) = s6*(T int ). From ( 1201) we obtain 



int 



6t + 6st + £ 3 t 2 + 0(| S ,t| 3 ) 
6r + ^r + e 3 -r 2 + 0(| S ,r| 3 ) 



0(2r) - s9{2t) < 
0(2r) - s9{2t) > 



(21) 
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where 



ad 1 G(e 1 ) 



-Wi sin(gi)G(gi) 
(a9 x G{9 x ) -sin(^)) 2 ' 



|6 sin 3 



(a0iG(0i) - sin(0i)) 3 ' 

P. = 1 fcCfg,) Sin2(gl) ~ 3Qgl g^lMglj + MiGjOi) 2 
3 2 (a9 1 G(9 1 ) -sin(6>i)) 2 

Note that the denominators of the £j share the same roots as (I14p when s = 0. As 
discussed in §2.3[ a root of (JHJ) corresponds to a boundary of a sliding region of the 
system in the absence of delay. 

From the expression for £i it follows that if a > g^gfgA , i-e. the control is sufficiently 
strong, and s and r are sufficiently small, then T indeed intersects S x at a time t = 
Ti n t > t. If G(6) = cos(d), this is equivalent to requiring a > 1 and 0q < 6* os . If 
G(6*) = 1, this condition is satisfied if a > 1 or 9\ < 9 < | . 

AH is derived by evaluating H(9(t), 4>(t)) at i = r and £ = T int + r, and taking the 
difference. The above expressions lead to 

AH-S ClST + C2T ' + &** T + C4Sr2 + Csr3 + T|4) ' Tint " 2T ^ 

I Cisr + C 2 r 2 + C 3 ^ 2 r + <>T 2 + C 5 r 3 + 0(\s, r| 4 ) , T int > 2r ' 
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where 

-a 2 0?G(0i) 2 



a9 1 G(9 1 ) - sin(0i; 



C2 = _ a ^ G{e ^)-h^K 



ae x G{0i) - sin(0i) ' 



6 _ ^j^H*^, . 

a9 1 G{9{) -sin(^i) X 

sin 3 ^) - ^sin^QG^) + 2a 2 g 2 sin^QG^) 2 - 
(,4 = — iao^G^ij -3 , 

f a0iG(0i) - sin(0i)J 

j. 1 ha ■ in ^ 2 sm 3 (fli ) - 6a6 1 sin 2 (gi)g(gi ) + 8a 2 fl 2 gMgi )G(9 l ) 2 - 3a 3 6fG(9 1 ) 3 
C5 = --awi sin(0iJG-(0iJ ^-3 ; 

b f a9 1 G{9 1 ) - sin(0i 

, , v9 sin 2 (0i) - |a0i sin(0i)G(%) + ±a 2 2 G(#i) 2 
U = 2ab9 2 l G{9 1 ) 2 — ^— V ' K ' x 4 2 — - — —— , 



a9 l G(9 l ) -sin(0i) 

1 n;a ^ 2 sin 3 (0i) + 7a0isin 2 (0i)G(ei) - 9a 2 2 sin^OG^i) 2 + 3a 3 0?G(#i) 3 
6 C 



C 5 = -abe x G{e x , 



a9 1 G{6i) -sm{6i\ 

From ( 1221) with Tj nt > 2r we may characterize zigzag dynamics described in the 
previous section (T int < 2r corresponds to relatively large values of a). T is a zigzag 
periodic orbit when AH = 0. The discontinuity-induced bifurcations of the previous 
section correspond to the creation of a zigzag periodic orbit at the origin. Therefore this 
bifurcation corresponds to AH = at an arbitrarily small value of 6\. Expanding AH in 
terms of 9\ allows us to determine the nature of the discontinuity-induced bifurcations: 

/ a 2 _ a 2 (a — 2) 2 ab(a — 2) 2 ab(a 2 — 3a + 4) 2 

V a-l 2(a-l) (a-1) 2 2(a - l) 2 

+ «fr(3a> - jg + 7a + 1) ^ + Q(|>[ r|4) j g; + 0(<?) (23) 

Note, (123]) is valid for both G(9) = cos(9) and G(9) = 1. The lowest order term (i.e. the 
9\ term) is zero when 

2 2(2 + 8a- 15a 2 + 6a 3 ) L 2 3 . , 0/l . 

-s i ^—^bs 2 + 0(s 3 ) . (24) 



a- 2 3a(a-l)(a-2) 3 

An omission of 0(s 3 ) terms in f [24"l) gives an approximation to the occurrence of the 
discontinuity-induced bifurcations and is shown in Figs. 0] and [5] (the dashed curves). 
As shown in these figures the approximation agrees well with the numerical results. The 
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approximation was obtained from series expansions in t and r and for this reason fits 
the numerics less precisely for values of a near 1, because here the period of the orbit is 
relatively large, and similarly worsens with increasing r. 

Since the second lowest order term in (1231) is of order two more than the previous term 
in 9, whenever this term is nonzero at a discontinuity-induced bifurcation the bifurcating 
zigzag periodic orbit grows at a rate proportional to the square-root of a non-degenerate 
change in the control parameters, a and b. This conclusion agrees with the bifurcation 
diagrams of Figs. H] and El Furthermore the criticality of the discontinuity-induced 
bifurcation is determined by the sign of the 9f term. Unlike the quadratic term, this 
term is dependent upon G"(0) and hence differs for G(9) = cos(#) and G(9) = 1. When 
G{6) = cos(9), the 9\ term vanishes when 

3 a _ 5 hi- 243a 6 + 1674a 5 - 4491a 4 + 5862a 3 - 3611a 2 + 642a + 175) „ , o N 
r s + S-7 TTT^ o 7vi L s 2 +0(s d ) 



3a 2 - 8a + 6 6a(a - 1) (3a 2 - 8a + 6) 3 

(25) 

Omitting 0(s 3 ) terms, fT23|) intersects the discontinuity- induced bifurcation curve (T2"4"l) 
at the point in Fig. [H-A indicated by a circle. Similarly when G(9) = 1, the 9\ term 
vanishes when 

t = _2 s _ 2^-2)(3.- 1) s2 + 
a 6a^{a — 1) 

The dashed curves of Figs. H] and |5] that approximate saddle- node bifurcations of the 
zigzag periodic orbits were computed numerically using quadratic and cubic terms of 

fl22D. 



3.5 Homoclinic bifurcations 

Here we derive, to lowest order in s and r, the curve of homoclinic bifurcations shown in 
Fig. HI At such a homoclinic bifurcation there exists the homoclinic connection shown 
in Fig. [71 This connection does not exist for relatively large values of r because in that 
case trajectories tend to intersect W s and undergo spiral motion. 

To approximate this connection we begin by deriving approximations to the stable 
and unstable manifolds of [9* os , 0). The eigenvalues and eigenvectors of the linearization 
about (#cos> 0) when r = are obtained by elementary calculations and give: 



focosf 



Ar 



± 



'b 2 cos 2 ( 



4 



+ 



29: 



sin(20* 



29* cosf 



(27) 
(28) 



where v ± is a vector in the (9, 0)-plane. The section of the orbit between the equilibrium, 
(#* os , 0), and the switching point, (92, ^2), see Fig. [JJ coincides with the unstable manifold 
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of (0* os , 0). From ( p7|) and (|28|) . this curve may be written as 

m = \+(9 - r m ) + o(\e- e* cos , r\ 2 ) . (29) 

We denote the intersection of ( J29l with Si (<fi = s9), by (6>i nt , s6> int ) and from (|29l) 
obtain 

0int = 0* cos (l + + 0(, 2 ) . (30) 

Near Ei the time-delay has more effect than near (#* os , 0). Nevertheless, from the same 
series expansion methods as the previous section, it follows that 

2 = ^nt + 0(\s,r\ 2 ) , (31) 

(i.e. the difference between 62 and 8 int is negligible for this calculation). From a series 
expansion in s and t of the solution to the OFF system between the two switching points 
(#2,^2) and (6x,<j>i) we obtain, using (13T|) . 

fli = flint + 0(|s,r| 2 ) , (32) 

0i = #intS + sin(6U)r + 0(kr| 2 ) • (33) 

Lastly, the section of the orbit between 0i) and (#* os , 0) is the stable manifold of the 
equilibrium so given by 

m = \^(6 - e: os ) + o(\e - e: os \ 2 ) . (34) 

The homoclinic connection exists when the point (#1, 0i), given by ( 13"i2"|) and f l3"3"|) . satisfies 
Using also (13"U|) . we arrive at (after some manipulation): 

' / 2 T + T+) s + °( s2 ) ■ ( 35 ) 



as a condition for the existence of the zigzag homoclinic orbit. The approximation 
obtained by dropping 0(s 2 ) terms in (|35l) is shown in Fig. H] and matches well to the 
numerical results. 




Figure 7: A sketch of a zigzag trajectory that is homoclinic to the equilibrium, (0* os , 0). 
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4 Simple and Complex Behaviour for Large Delay 



Dynamical behaviour exhibited by the system for small delay was described in the 
previous section. However we did not provide a complete description of dynamics in 
the case that a is slightly larger than 1. In this case complex behaviour may occur 
that appears to persist for larger r. In section §4.11 we demonstrate that solutions of 
the system may exhibit distinct behaviours on different time-scales in a manner akin to 
bursting in neuron models. In section §4.21 we describe four bifurcations that indicate 
behaviour near the origin when the delay time is large. 

4.1 Bursting-like dynamics 

Fig. |8]-A illustrates dynamics of (HJ)-([6]) when r = 0.3 and s = —0.1 for different values of 
a. The majority of the dynamics indicated by this plot match the predictions of Fig. H] 
for behaviour for small r. Specifically there are stable and unstable zigzag periodic 
orbits that grow in size proportional to the square root of change in a and that collide 
and annihilate in a saddle-node bifurcation. However for a% < a < 03 (where a<i ~ 1.161 
and 03 ~ 1.208, see Fig. |8]-A), numerical simulations reveal a complicated attracting set. 
This set is shown in panels B and C of Fig. [8] for a = 1.18. 

To describe this set and its associated dynamics consider the forward orbit of a 
point on Si, near the origin. For the parameter values of Fig. [8]-B, the forward orbit 
zigzags away from the origin, for some time. Over the course of each zigzag, the orbit 
spends a time greater than r within the OFF region that decreases for each successive 
zigzag. The combination of slowly outward-moving motion (with 9 = 0(r)), and rapid 
zigzag motion (of frequency 0(~)) continues until 9 ~ 0.32 at which the orbit enters 
and exits the OFF region in a time less than r. Let t s hort < t denote this time and let 
(#short, s^short) denote the point on Si at which the orbit exits the OFF region. Unlike 
what has been typical throughout this paper, beyond (0 s hortj s^short) the orbit is governed 
by the ON system for the next r — t s h or t time after which it is governed by the OFF 
system for the next t s hort time. Numerically we observe that because this part of the 
orbit follows the OFF system for less time than it would on a typical zigzag oscillation, 
it does not attain such a large Rvalue before the control is reapplied. Consequently the 
orbit may to become trapped in the ON region, albeit only slightly above Si, for some 
time. This is seen in Fig. [8]-C. During this time the orbit rapidly approaches the unstable 
manifold of the saddle equilibrium, (0* OB ,O). The unstable manifold tends to the origin 
but for the parameter values of Fig. [8]-B, the manifold intersects Si at 9 ~ 0.236. This 
#-value is sensitive to the choice of the value of a, as visible in Fig. [8|A. 

Therefore, after a zigzag oscillation involving a time in the OFF region that is less 
than the delay time of the system, r, the forward orbit endures an excursion back to Si 
(at 9 ~ 0.236 in Fig. [8]-B) followed then by outward- moving zigzag motion and continued 
repetition of this procedure. For a = 1.18 the attracting set is periodic but for different 
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values of a with a 2 < a < a^, the attracting set may be chaotic. 

As the value of a is decreased from a = 1.18, the range of values over which the 
attracting set exists, increases, until at a = a 2 the unstable manifold of (0* os ,0) no 
longer intersects £]_. For a < a 2 trajectories following the manifold limit directly to the 
origin. It is interesting to note that for a\ < a < a 2 (a% ~ 1.135), the forward orbit of 
a point on Si arbitrarily close to the origin zigzags slowly away from the origin until 
at a #-value of order 1 the orbit becomes trapped in the ON region and approaches the 
origin. Consequently over this range of a-values the origin is not Lyapunov stable, and 
hence not asymptotically stable jH], but appears to be quasi-asymptotically stable in 
that all points in a neighbourhood of the origin tend to origin, eventually. 

The complicated attracting sets are born out of the stable zigzag periodic orbit in 
a bifurcation at a = 03, Fig. [8]- A. At the bifurcation the amount of time spent by the 
periodic orbit in the OFF region is exactly the delay time, r. This type of discontinuity- 
induced bifurcation was analyzed in [21] (see also for a general time-delayed, 
piecewise-smooth system comprised of ordinary differential equations on each side of 
the switching manifold. In that paper it was proved that in a neighbourhood of the 
bifurcation a Poincare map is generically piecewise-smooth continuous, and to lowest 
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Figure 8: Panel A is a bifurcation diagram of (SJ)- (J6j) when r = 0.3, s = —0.1, 6 = 2 
and G{9) = cos(0). The abbreviations and line styles are explained in the caption of 
Fig. 0]-B. Panel B shows a partial time series of an orbit with transients decayed when 
a = 1.18 (corresponding to the dash-dot line in panel A). Panel C shows this orbit in 
the (9, 0)-plane. For values of a near 1.18, in panel A we have indicated ^-values at 
which the orbit crosses Si and then immediately enters the neighbouring ON region. 
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order piecewise-linear. However, for the system studied in this paper phase space is 
infinite-dimensional so it is not clear how to define a Poincare section that captures 
all oscillatory motions local to the bifurcation. We leave for future work a thorough 
investigation of the discontinuity-induced bifurcation characterized by a time spent in 
the OFF region exactly equal to the delay time of the system. 



4.2 Four fundamental bifurcations for dynamics near the origin 

It is more difficult to classify dynamics of the system when the delay time, r, is large. 
For large r, spiral dynamics (described in §3.ip may dominate. Spiral dynamics cannot 
be analyzed by the asymptotic methods of §3.41 because the time taken for a single 
spiral is order 1. Furthermore, spiral periodic orbits may undergo symmetry breaking 
bifurcations followed by period-doubling cascades to complex attractors. We leave a 
more complete analysis of these bifurcations for future work. In this section we analyze 
the stability of the origin. 

The origin is a non-differentiable point of (jl])-© and so does not have associated 
eigenvalues that determine stability. The stability of such a point in a piecewise-smooth 
ODE system is understood in two dimensions but is yet to be completely solved in 
higher dimensions [4"2"| |4"3"| HU US] . The presence of time-delay only adds complexity; for 
this reason we rely on numerical simulations. Since the sole interest here is on dynamics 
local to the origin, it is sufficient to analyze the linearization of (jl]): 

e = </>, 

. (36) 
(j) = 6 + F , 

which is valid for both G(6) = cos((9) and G{6) = 1. Due to the scale invariance of the 
spatial coordinates it suffices to look at the forward orbits of just one point on Si, say 
(1, s), and one point on S 2 , say (0, 1). 

With the goal of determining the dynamics for all combinations of the control pa- 
rameters, we perform numerical integration to study these two forward orbits. For each, 
we identify, and numerically continue, two key bifurcations. First, the forward orbit of 
(l,s) may return to this point upon one zigzag. In the (a, &)-plane this occurs along 
the solid curve in Fig. [9j for which r = 0.5 and s = 0. We have found that different 
values of r and s produce qualitatively similar pictures. Locally, for values of a and b 
to the left of this curve, the forward orbit of (1, s) zigzags away from the origin; to the 
right of the curve the orbit zigzags into the origin. For the nonlinear system, this is 
the discontinuity- induced bifurcation identified in §3.31 at which two symmetric zigzag 
periodic orbits are created. Second, the forward orbit of (1, s) may become coincident 
to W s (the stable manifold of the origin for the OFF system, Fig. [I]) . This occurs along 
the dashed curve in Fig. For values of a and b to the left of this curve the forward 
orbit zigzags, to the right of the curve it spirals. 
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For the initial point (0, 1) there are two bifurcations analogous to those just discussed. 
Along the dash-dot curve in Fig. [9] the forward orbit of (1, s) returns to (1, s) after one 
spiral, along the dotted curve the orbit falls onto W s and limits upon the origin without 
again crossing either of the switching manifolds. The dash-dot curve is a bifurcation 
of (jl]) at which a symmetric spiral periodic orbit is born in a manner akin to a Hopf 
bifurcation. 

Note that one may choose the control parameters such that zigzag orbits approach 
the origin and spiral orbits head away from the origin (as in Fig. [1]) and vice-versa. We 
have not been able to identify an intersection between the dashed and dotted curves 
of Fig. ED for any r and s, nor have been able to show that such an intersection cannot 




Figure 9: A bifurcation set of the linearization, ()36l) . with (EJ)- (J6j) , r = 0.5 and s = 0. 
The (a, 6)-plane has been partitioned according to the fate of forward evolution of the 
points (1, s) and (0, 1). By zigzag in [zigzag out] we mean that the forward orbit of (1, s) 
zigzags into [away from] the origin. Similarly by spiral in [spiral out] we mean that the 
forward orbit of (0, 1) spirals into [away from] the origin. By zigzag to spiral we mean 
that the forward orbit of (l,s) undergoes spiral motion and behaves like the forward 
orbit of (0, 1), and vice- versa for spiral to zigzag. For instance when (a, b) = (3, 3) both 
forward orbits spiral away from the origin; when (a, b) = (1.5, 3.5), the forward orbit of 
(1, s) zigzags in to the origin and the forward orbit of (0, 1) spirals away from the origin, 
as in Fig. [IJ For comparison we have also indicated, by the gray D-shaped region, the 
region where the origin is a stable equilibrium of the ON system. Within this region the 
countable set of eigenvalues associated with the equilibrium all have negative real part 

OUT]. 
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occur. Such an intersection could permit for the existence of an orbit that repeatedly 
switches between zigzag and spiral motion. 

5 Dynamics with the Switching Condition (JTj) 

In this section we analyze the system ([I])-© with the alternative switching rule dZj), 
Fig. [TUJ With this rule the control is removed when the time-delayed position is near 
vertical and is motivated from observations of human balancing tasks and postural sway 
[T9] , Switching control off near the origin may lessen "over-control" but eliminates the 
possibility of a stable vertical position. Previous investigations have used equations of 
motion that are linear in 6 [191 EH HH]. For typical practical applications this is justified 
because the relevant range of 6 values is sufficiently small. We have chosen not to 
linearize in 6 since our asymptotic methods do not rely on it and so that we may study 
the influence of additional equilibria. 

The system fll])-© with ([7]) is infinite-dimensional, so, as above, we look only at the 
forward orbits of points on switching manifolds. We believe that an analysis of the fate 
of these orbits provides a good understanding of the important properties of the system 
for a wide range of parameter values. Orbits that cross the OFF region in a time less 
than t occur readily if r is large relative to the width of the OFF region. In this case 
typical stable dynamics are oscillations that involve a relatively large range of ^-values 
which spend only a small fraction of time in the OFF region and are therefore often 
similar to dynamics of simply the ON system. 
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Figure 10: The (9, 0)-plane for ®-© with © and a = 0.3, r = 0.5, (a, b) = (2.5,4) 
and G(9) = cos(#). A trajectory approaches an asymmetric periodic orbit. W s and W u 
are the stable and unstable manifolds of the origin for the OFF system, respectively. 
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Figure 11: Bifurcation diagrams and representative sketches of dynamics in the (9,4>)- 
plane of the system ®-© with (J7J) for r = 0.1, a = 0.3, b = 0.5 and G(6) = cos(0) 
in panel A and G{6) = 1 in panel B. HC - homoclinic bifurcation; BEB - boundary 
equilibrium bifurcation. Due to the symmetry of the system, dynamics for 9 < are 
identical to those for 9 > and for this reason are not shown. Between the homoclinic 
bifurcations the double curves indicate the maximum and minimum ^-values of a sta- 
ble periodic orbit. To the right of the right-most homoclinic bifurcation there exists 
one symmetric periodic orbit; its minimum value is not visible in the bifurcation dia- 
grams. Solid [dashed] curves correspond to stable [unstable] equilibria. The dotted lines 
represent the switching manifold, 9 = a. 
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For ease of explanation we discuss dynamics only for 9 > 0; by symmetry identical 
dynamics occur for 9 < 0. On the switching manifold S3 t^j, 9 = <j>, thus trajectories 
that cross £3 at, say, (er, <p ), next enter the neighbouring ON region if </> > and next 
enter the OFF region if <j) < 0. For this reason we study the forward orbits of points 
(er, O ) with 0o > in view of our earlier discussion regarding initial conditions, §2.21 

Fig. [Til shows bifurcation diagrams when r = 0.1, a = 0.3 and b = 0.5. Note that 
the value of a used in this illustration is significantly larger than values suitable for 
traditional human balancing tasks so may be more meaningful in regards to mechanical 
applications. However, the purpose of Fig. [TTJis to highlight the basic bifurcation struc- 
ture of the model. We have identified qualitatively similar dynamics over a wide range 
of parameter values, in particular with smaller values of er and r. The most parame- 
ter sensitive component of the bifurcation structure is the value of a of the right-most 
homoclinic bifurcation which decreases with an increase in r. 

The equilibrium of the ON system, (#* os ,0) or (91,0), is admissible if and only if 
it lies to right of £3. Upon variation of a, the equilibrium collides with the switching 
manifold at the point (er, 0). (This is referred to as a boundary equilibrium bifurcation 
|39|.) For values of a between the two homoclinic bifurcations, Fig. [TU there exists a 
stable periodic orbit encircling (er, 0) (and a symmetric orbit encircling (— er, 0)). As the 
value of a is decreased the periodic orbit is destroyed in a homoclinic bifurcation with 
(9* os , 0) or (91, 0) (depending on the function G{9)). In the case G(9) = 1 this bifurcation 
is coincident with the boundary equilibrium bifurcation. As the value of a is increased 
the two symmetric periodic orbits connect at the origin in a homoclinic bifurcation 
beyond which there exists one stable symmetric periodic orbit encircling the origin. 
With a further increase in a, or an increase in r, numerically we have observed that 
this periodic orbit undergoes symmetry breaking and period doubling to an aperiodic 
attractor in a fashion similar to the system with continuous control [12]. 

For the remainder of this section we perform an asymptotic analysis, along the same 
lines as in §3.4l for fl6]), to analytically derive the small amplitude periodic orbit encircling 
(er, 0) for small r, and determine the rate at which the periodic orbit grows in size with 

T. 

For small O > 0, let T be the forward orbit of (a, (f> ) at t = 0, for (TjJ-flS]) with 
([7]) and, recalling the discussion in §2.2[ assume that T lies in the OFF region for all 
t G [— r, 0]. Switching occurs within the ON region at t = r and, assuming T re-enters 
the OFF region at some later time T int , a second switching occurs at t — T int + r before 
T exits the OFF region. We let and (#21^2) denote the respective switching 

points, Fig. [T^J Our goal is to calculate the change in the Hamiltonian ( II Op between the 
two switching points in order to identify periodic orbits. Let 

= 00 , H 2 = H(9 2 ,cp 2 ) . 
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r is periodic if AH = H2 — Hi = 0. Since Hi = H(a, 4>o), we have 

1 



Hi = 



cosier 



(37) 



As in §3.41 we obtain a useful description of T by substituting a series representation 
of 8(t) expanded in 0o and t into the equations of motion and solving for the unknown 
coefficients: 



6(t) 



( a + <pQt + \ sin(cr)i 2 + O{\<f> , t| 4 ) , t € [0, r] 

(a + o r + I sin((r)r 2 ) + (0 O + sin(cr)r)(t - r) 

+ (a 3 + f0o)(t-r) 2 + ^(t-r) 3 + O(|0o, r,t| 4 ) , t G [r, 2r] 

(a + 4> t + i sin((r)r 2 + air 3 ) + (<p + sin(cr)r + d 2 T 2 )(t - r) 

+ (a 3 + + a S r)(t - r) 2 + a 6 (t - r) 3 + O(|0 O , r, t| 4 ) , t G [2r, T int + r] 



(38) 

where the a's, listed in §3.4[ are evaluated at 6 = a. To determine Tj nt from 9(T int ) = a 
it is necessary to consider 4> = 0(ra) and write 

Tint = + + X3T + O(|0o, T^| 3 ) . (39) 

Analogous to §3.4[ the unknown coefficients are calculated by substituting fl39|) into 
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Figure 12: A sketch of the forward orbit of a point on E 3 for the system flU)-© with 
the switching rule (jTJ). 
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©8]): 

Xi 

X2 
X3 



aaG(a) — sin(cr) ' 
-I' 



}bG(a) 



{aaG{cx) - sin(a)) 2 ' 

2aaG(a) 
aaG(a) — sin(a") 



assuming T int > 2r. Notice we must have sin(cr) — aaG{a) < because we require 
T int > 0. An evaluation of (fTU|) at t = T int + r using (157)) . (155)) and (15il produces 

-bG(cr} 

AH = 2aaG(a)<j) T -P j ^ + O(|0 O ,^| 4 ) . (40) 

aaG{a) — sm(crj 

Consequently AH = when 



T 



3aa(aaG(a) — sin(er)) 



<Pl + 0{4>l), (41) 



which is consistent with our assumption that 0o = 0(t*). 

In summary, when AH = 0, T is a small amplitude periodic orbit encircling ((7,0), 
see for instance Fig. [TTJ As the value of r is increased from zero, we deduce from 
(I4ip this periodic orbit grows out of the point (er, 0) with an amplitude asymptotically 

proportional to r^. Furthermore, at this periodic orbit = aaG ^y}^l^ + ^(^o)' 
which is negative-valued because aaG(a) — sin(cr) > 0. Consequently the periodic orbit 
is stable matching the numerical results of Fig. [TTJ 



6 Discussion 

In this paper we have identified bifurcations and dynamics in a prototypical balancing 
model describing planar motion of an inverted pendulum with control that is qualita- 
tively affected by the combination of time-delay, discontinuity in the control, and nonlin- 
earity. Time-delay is fundamental to a variety of balancing problems. It represents the 
time-lag between when variables are measured and corrective forces are applied, which 
in human balancing tasks typically represents neural transmission time. Switching in 
the method of control has been proposed to reflect observations of intermittent mus- 
cle movements, to procure simplicity in mechanical systems, or to provide a stabilizing 
mechanism particularly when the time-delay is long. Finally, terms in the equations of 
motion that are nonlinear in the angle of displacement from vertical, 9, are important 
when the value of 9 is not restricted to small values. 
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Previous work uses mathematical methods to analyze systems that lack any one of 
these three features. The bifurcation theory and methods of piecewise-smooth systems 
[3D] apply to systems without time-delay. Centre manifold reductions may be applied 
to models that lack a switching condition and are smooth [To"! [T7] . Systems that are 
spatially scale-invariant have been considered in the context of balancing [22J and in 
general (121 H31 HU US]; numerical simulations are often essential in this case. Even 
though (jH)-© with either (jSJ) or ([7j) exhibits all three of the above features we have 
been able to obtain some formal results. One simplifying aspect is that the system does 
not switch between two DDEs, as in for instance [2H EH], but rather switches between a 
DDE (the ON system) and an ODE (the OFF system). As a result, whenever an orbit 
spends a continuous length of time equal to or greater than the delay time, r, governed 
purely by the OFF system, its future evolution becomes independent of its location at 
any earlier time. Consequently initial conditions of the system may be thought of as 
points in the (9, 0)-plane. More specifically, since the dynamics of the OFF system is 
lucid, for initial conditions we use points on the boundaries of the ON/OFF regions at 
which the vector field of the OFF system points into the neighbouring ON region, §2.21 

The presence of time-delay in the switching rules induces what we have referred to 
as zigzag motion. This motion is characterized by a rapid on/off switching of the control 
and corresponds to a restriction of the pendulum to one side of the vertical position. 
Since zigzag motion occurs on an 0(r) time-scale, it succumbs to the asymptotic ap- 
proach, based on piecewise Taylor series in t and r, for example (1201) . We performed the 
asymptotic analysis for the particular state-dependent switching rule, (J6]), where we ex- 
panded also in the switching parameter s. When — s > r, and r is not too large, zigzag 
motion of the pendulum approaches the vertical position on a relatively long time-scale, 
§3.21 This manner of stabilizing the vertical position is not possible without a switching 
rule like <Q. 

Nonlinearity in 9 in (jl])-© with ([6]) permits non-equilibrium, asymptotically stable 
invariant sets. Using the series expansions mentioned above, we have been able to 
identify periodic orbits of period 0(r), and derive equations in terms of the system 
parameters that correspond to bifurcations of these periodic orbits. For instance zigzag 
periodic orbits bifurcate from the vertical position in symmetric discontinuity-induced 
bifurcations. We have analyzed the stability of these periodic orbits both numerically, 
§3.3[ and through asymptotic expansions, §3.41 A homoclinic bifurcation was identified 
for small r and s in a similar fashion, §3.51 We also described a complicated bursting-like 
attractor in §4.11 

For relatively large values of r the model predicts the pendulum to typically prefer 
oscillations about the vertical position. We have referred to such motion as spiral motion 
due to the nature of corresponding trajectories in the (9, </>)-plane. In §4.21 we have 
investigated the spiral motion numerically in the context of the stability of the vertical 
position. Since the spiral behaviour operates on long time-scales it cannot be analyzed 
by the asymptotic approach described above. In particular we found that the vertical 
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position may be semi-stable in that for a fixed choice of control parameters, spiral motion 
may approach the vertical position whereas zigzag motion heads away from this position, 
or vice-versa. It is interesting to note that dynamics local to the vertical position are 
explained by a global analysis of the piecewise-linear system, ( I36p . This circumstance of 
global dynamics governing local behaviour is a common occurrence in piecewise-smooth 
systems, see for instance |46j . 

In §5]we studied the model with the switching rule (J7J) that turns off the control when 
the controller interprets the magnitude of 9 to be less than a threshold value, a. In this 
setup the vertical position is always unstable. Using the same asymptotic methods we 
showed that as the value of r is increased from zero a stable periodic orbit emanates 
from (9, (f>) = (a, 0) with an amplitude asymptotically proportional to ra. This periodic 
orbit also corresponds to small periodic fluctuations of the pendulum on one side of 
the vertical position due to an intermittent application of the control. By symmetry 
there exists an identical stable periodic orbit on the other side of the vertical position. 
As the value of r is increased, typically the two periodic orbits collide in a homoclinic 
bifurcation with the vertical position beyond which there exists one symmetric stable 
periodic orbit corresponding to oscillations about the vertical position. With a further 
increase in r dynamics exhibited by the system are similar to the case where the control 
is constantly applied. 

A future project is that of an investigation into the effect of noise in models of the 
type studied here. Some steps in this direction have already been achieved [22j [301 HZ]- 
Noise may result from discrepancies in measurements of controller or from fluctuations 
in muscle response and may induce a flip-flop motion between coexisting stable solutions 
or possibly have a stabilizing effect [HI HE] . 
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